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Abstract 

We study the approach to jamming in hard-sphere packings, and, in particular, the pair corre- 
lation function g2{r) around contact, both theoretically and computationally. Our computational 
data unambiguously separates the narrowing delta-function contribution to g2 due to emerging 
interparticle contacts from the background contribution due to near contacts. The data also shows 
with unprecedented accuracy that disordered hard-sphere packings are strictly isostatic, i.e., the 
number of exact contacts in the jamming limit is exactly equal to the number of degrees of freedom, 
once rattlers are removed. For such isostatic packings, we derive a theoretical connection between 
the probability distribution of interparticle forces Pf{f), which we measure computationally, and 
the contact contribution to g2- We verify this relation for computationally-generated isostatic 
packings that are representative of the maximally jammed random state. We clearly observe a 
maximum in Pf and a nonzero probability of zero force, shedding light on long-standing questions 
in the granular-media literature. We computationally observe an unusual power-law divergence 
in the near-contact contribution to g2, persistent even in the jamming limit, with exponent —0.4 
clearly distinguishable from previously proposed inverse square root divergence. Additionally, we 
present high-quality numerical data on the two discontinuities in the split-second peak of g2, and 
use a shared-neighbor analysis of the graph representing the contact network to study the local 
particle clusters responsible for the peculiar features. Finally, we present the first computational 
data on the contact-contribution to 52 for vacancy-diluted FCC crystal packings and also investi- 
gate partially crystallized packings along the transition from maximally disordered to fully ordered 
packings. Unlike previous studies, we find that ordering has a significant impact on the shape of 
Pf for small forces. 



* Electronic address: |torquato@electron.princeton.ei 
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I. INTRODUCTION 



Jamming in hard-sphere packings has been studied intensely in the past years. In this 
paper, we investigate the pair correlation function (72 (r) of the classical three-dimensional 
hard-sphere system near a jamming point for both disordered (amorphous, often called 
random) as well as ordered (crystal) jammed packings. The basic approach follows that of 
Ref. Q], developed further for crystal packings of rods, disks and spheres in Ref. |2]. We 
focus on finite sphere packings that are almost collectively jammed ^|4|, in the sense that 
the confirmation pent is trapped .n a very small region of eonflguration spaee around the 
point representing the jammed ideal packing Difficulties with extending the results to 
infinite packings will be discussed in what follows. In the ideal jammed packing particle 
contacts necessary to ensure jamming are exact, and the particles cannot at all displace, 
even via collective motions. Such ideal jammed (or rigid) packings have long been the 
subject of mathematical inquiry jsj; however, they are not really attainable in numerical 
simulations where produced packings invariably have some interparticle gaps (even taking 
into account the unavoidable roundoff errors) . It is therefore instructive to better understand 
the approach to this ideal jammed state computationally and theoretically, which is the 
primary objective of this paper. 

We choose as our main tool of exploration the shape of the venerable orient at ionally- 
averaged pair correlation function g2{r) around contact. This is because this function is 
a simple yet powerful encoding of the distribution of interparticle gaps. In the jamming 
limit, it consists of a delta function due to particle contacts and a background part due to 
particles not in contact. As the jamming limit is approached, it is expected that the delta- 
function contribution will become more localized around contact. We derive the first exact 
theoretical model for this narrowing for isostatic packings (defined below), connecting (72 to 
the probability distribution of interparticle forces Pf, and verify the relation numerically. 
In this work, we present computational data with unprecedented proximity to the jamming 
limit, for the first time clearly separating the narrowing delta-function contribution from the 
apparently persistent diverging background contribution. The data show that our disordered 
packings have an exactly isostatic contact network in the jamming limit, but with an unusual 
multitude of nearly closed contacts. We study the properties of the contact network and find, 
contrary to previous studies, no traces of polytetrahedral packing, but rather a complex local 



3 



geometry, indicating that the geometric frustration due to the constraints of global jamming 
on the local geometry is nontrivial. 



II. THEORY 



A packing of hard spheres of diameter D in c?-dimensional Euclidian space is charac- 
terized by the (A^(i)-dimensional configuration vector of centroid positions R = (ri, . . . , r^r). 
Here we fix the center of mass of the packing (with periodic boundary conditions), so that 
in fact the configuration space is of dimension (A^ — l)d. However, we will usually neglect 
order unity terms compared to A^. The boundary conditions imposed determine the volume 
of the enclosing "container" V and the packing (covering) fraction, or density, 0. 

A jammed packing is one in which the particle positions are fixed by the impenetrability 
constraints and boundary conditions 0, Q]. In particular, a packing is locally jammed if 
no particle in the system can be translated while fixing the positions of all other particles; 
collectively jammed if no subset of particles can simultaneously be continuously displaced 
so that its members move out of contact with one another and with the remainder set; and 
strictly jammed if it is coUectivelly jammed and all globally uniform volume- nonincreasing 
deformations of the system boundary are disallowed by the impenetrability constraints. 
Assume that a configuration Rj represents a collectively jammed ideal packing :4] with 
packing fraction (pj, where there are M interparticle contacts. Next, decrease the density 
slightly by reducing the particle diameter by AD, 6 = AD/D ^ 1, so that the packing 
fraction is lowered to = 0j (1 — 6)'^. In this paper we restrict ourselves to an analysis 
which is first order in the jamming gap 5, ~ 0j (1 — d6), and focus on three-dimensional 
packings, d = 3. 



A. Jamming: Configurational Trapping 

It can be shown that there is a sufficiently small 5 that does not destroy the jamming 
property, in the sense that the configuration point R = Rj + AR remains trapped in a 
small neighborhood JTar around Rj In fact, for sufficiently small S, it can be shown 
that asymptotically the set of displacements that are accessible to the packing approaches 
a convex limiting polytope (a closed polyhedron in arbitrary dimension) Par ^ iTar 



4 




Figure 1: The polytope of allowed displacements "Par for a locally jammed disk (light shade) 
trapped among three (left) or six (right, as in the triangular lattice) fixed disks (dark shade). 
The exclusion disks (dashed lines) of diameter twice the disk diameter are drawn around each of 
the fixed disks, along with their tangents (solid lines) and the polytope T^ar they bound (dark). 
For the isostatic case on the left this polytope is a triangle (a simplex in two dimensions), and a 
hexagon for the hyperstatic case on the right. 

This polytope is determined from the linearized impenetrability equations 

A^AR < Al, (1) 

where A is the (dimensionless) rigidity matrix^ of the packing, and Al is the set of inter- 
particle gaps [4j. In our case Al = ADe, where e is a vector of M elements all equal to one. 
We can therefore focus on the normalized polytope • A-^x < e, which can be scaled by 
a factor of 6D to obtain Par- Examples of such polytopes for a single disk are shown in 
Fig. A troublesome aspect, discussed in Ref. is that infinite packings can never be 
jammed in the above sense unless 6 = 0, due to the appearance of unjamming mechanisms 

^ This matrix combines geometrical information with the topological connectivity information contained in 
the node-arc incidence matrix of the graph representing the contact network of the packing. Namely, A 
has Nd rows, d rows for each particle, and M columns, one for each contact. The column corresponding 
to the contact between particles i and j is nonzero only in the rows corresponding to the two particles, 
and contains the unit surface normal vector at the point of contact 0| . 
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involving collective density fluctuations. Nevertheless, computational studies indicate that 
macroscopic properties derived using this polytope-based approach do not depend on A^, 
even as — ^ oo. 

The polytope Px is necessarily bounded for a jammed conflguration, which implies that A 
is of full rank [J], and that the number of faces bounding Vy^, i.e., the number of interparticle 
contacts M, is at least one larger than the dimensionality dcs of the conflguration space, ^ 
M > dcs + ^- For collective jamming the boundary conditions are fixed and with periodic 
boundary conditions there are d trivial translational degrees of freedom, so dcs = {N — l)d. 
If hard wall boundary conditions are employed then dcs = Nd and one should also count 



contacts with the hard walls among the M constraints. For strict jamming ^ the boundary 
is also allowed to deform and this introduces additional degrees of freedom. For example, 
with periodic boundary conditions a symmetric non-expansive macroscopic strain tensor is 
added to the conflguration parameters, giving dcs = {N — l)d + d{d — l)/2 + {d — 1) degrees 
of freedom^. Isostatic packings are jammed packings which have the minimal number of 
contacts, namely, for collective jamming 

f 2A^ - 1 for d = 2 
M=l (2) 

I 3iV - 2 for = 3 

and for strict jamming 

i 2N + lioT d = 2 
M = { , (3) 

3A^ + 3 for = 3 

with periodic boundary conditions. Packings having more contacts than necessary are hyper- 
static, and packings having less contacts are hypostatic (for sphere packings these cannot be 
jammed in the above sense). For the trivial example of local jamming and = 1, all parti- 
cles but one are frozen in place and the free particle must have at least d + 1 contacts. Figure 
shows the polytope Par for a locally jammed disk, for both an isostatic and a hyperstatic 
case. In this work, we focus on collectively jammed packings, since strictly jammed packings 



^ The additional +1 comes because we are considering inequality constraints, rather than equalities. One 

can also think of this extra degree of freedom as representing the density, i.e., the size of the particles. 

For example, looking at the left panel of Fig. ^we see that at least 3 linear inequalities are necessary to 

bound a polytope in 2 dimensions. 
^ Here d{d — l)/2 gives the number of off-diagonal strain components, and d~ 1 comes from the number of 

diagonal components {d) whose sum is constrained to be non-positive (—1). 
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are hard to produce with existing algorithms. For sufficiently large disordered systems, the 



0. 



differences between collective and strict jamming are expected to be insignificant 

We now consider adding thermal kinetic energy to this nearly jammed hard-sphere pack- 
ing. While the system may not be ergodic and thus not in thermodynamic equilibrium, 
especially if considering disordered packings Q], one can still define a suitable macroscopic 
pressure by considering only time averages as the system executes tightly confined motion 
around the particular configuration Rj. In a sense, the configuration will explore the interior 
of "Par, and ergodicity is restored if one restricts the configurational space to "Par- For a 
finite packing, which is sufficiently close to its jamming point, the time-averaged properties 
will always be well-defined. Since the available (free) configuration volume scales in a pre- 
dictable way with the jamming gap, I^arI = (SD)^'^ ["Pxl, one can show that the reduced 
pressure is asymptotically given by the free- volume equation of state 

= ^ = -= ^ (A) 

^ NkT 6 (1-0/0J)' 

Relation Q is remarkable, since it enables one to accurately determine the true jamming 
density of a given packing even if the actual jamming point has not yet been reached, just 
by measuring the pressure. We later numerically verify the validity of Eq. Q in the vicinity 
of the jamming point. 

B. Jamming: Interparticle Forces 

As the particles travel around Rj and the configuration explores "Par? one can average 
the exchange of momentum between any two pairs of particles which share a contact in the 
jammed limit (i.e., whose contact forms a face of Vx), hereafter referred to as first neighbors, 
to obtain an average interparticle force (momentum transfer per unit time) |9| . The vector 
of coUisional forces f compares directly to the inter-grain force networks which have been 
the subject of intense experimental and theoretical study in the field of granular materials 
[10, luJ, ll3 ! even though we are to our knowledge the first ones to directly observe and 
measure them for true hard particles j^. These forces are in local equilibrium, 

Af = 0, (5) 
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where we take the forces to be nonnegative^, f > 0, and normahze them to have a unit 
average, / = e^f/M = 1, in the tradition of the granular media hterature. Our numerical 
investigations indicate that indeed the set of time-averaged collisional forces approach local 
equilibrium as the time horizon T of the averaging increases, in a inverse-power law man- 
ner, ||Af II ~ T^^. We can therefore obtain interparticle forces relatively accurately given 
sufficiently long molecular dynamics runs. While Eq. (jSJ will have a unique solution if and 
only if the contact network of the packing is isostatic, even for hyperstatic packings, such 
as the FCC packing, the equilibrium set of forces should be unique. In fact, one can prove 
that the force between two particles will be proportional to the surface area of the face of 
Px formed by the contact in question. 

It is interesting to observe that if one has an arbitrary point AR G Par? the interparticle 
gaps due to nonzero jamming gap will be Al ^ A"^AR — 6De, so that 

f^Al ^ (Af )^AR - MD6 = -MD6. (6) 

Eq. © enables one to determine how far from the jamming density a packing is without 
actually reaching the jamming point. This can be a useful alternative to using Eq. (0} when 
the hard-sphere pressure is not available, but interparticle forces are, such as, for example, 
with packings generated by algorithms using stiff "soft" spheres . 

As already pointed out, hypostatic packings cannot be jammed. However, it is possible 
for a hypostatic packing to be locally maximally dense, in the sense that no continuous 
motion of the particles can increase the density to first order. In other words, the particles 
must first move and unjam (which must be possible for a hypostatic sphere packing) before 
the density can increase. In particular, a packing of contacting particles for which a set 
of interparticle forces f in equilibrium exists, is locally maximally dense. In a sense, the 
interparticle forces resist further increase of the density. As we discuss later, our packing 
generation algorithm sometimes terminates with such packings since it tries to continually 
increase the density. 

This sign convention is in agreement with the granular media hterature, but opposite to our own preferred 
notation Ijl. 
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C. Pair Correlation Function Around Contact 



We now turn to the central subject of this work: The shape of the (orientation-averaged) 
pair correlation function (72 (r) for small jamming gaps. In particular, we will focus on 
interparticle distances r that are very close to D. We express (72(0 i'^ terms of the nonnegative 
interparticle gaps I = r — D. The polytope picture above says that only the M first-neighbor 
particle pairs will contribute to the shape of g2{l) right near contact, i.e., for gaps up to 
^max, where /max IS the largest distance from the centroid of Par to one of its faces. This 
contribution will become a delta function in the jamming limit. Particle pairs not in contact 
will not contribute to (?2(0 until gaps larger than the minimal further-neighbor gap IpN, 
and for now we will implicitly assume that Ipjy ^ /max, so that there is a well-defined 
delta-function region g2\l) = fi'2(/ <^ Wn)- This delta function region has previously been 
investigated theoretically for crystal packings, primarily j2j- In this work, we derive exact 
theoretical expressions for this region for isostatic packings, as well as numerically study 
vacancy-diluted FCC crystals and partially crystallized packings. 

1. Isostatic Packings 

We first focus on the probability distribution for observing an interparticle gap /, Pi(/), 
which is related to g2\l) via a simple normalization factor. The contribution P{1) from 
a specific contact is determined from the area S{1) of the cross section of with a plane 
)arallel to the face corresponding to the contact and at a distance / from the face, P{1) ~ S{1) 
^. The critical observation we make is that for an isostatic contact network, Px is a simplex 
and thus immediately we get S{1) ~ [{h — l)/Kf'^, where h is the height of the simplex 
corresponding to this particular face, h = M \Vx\ /S, S = 5(0). After normalization of P{1) 
and averaging over all interparticle contacts, we obtain that 



which shows that if we know the distribution Ph of heights for the simplex Px, or equivalently, 
the distribution of surface areas S of the faces of the polytope Ps{S), we would know Pi and 
thus g2^. 

Since the interparticle force / ~ S*, we see immediately that the distribution of face areas 
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is equivalent to the distribution of interparticle forces Pf{f), and in fact it is easy to derive 
that 



which gives in the hmit M oo 



m) 







lf=o 


MAD 



M 



Pf{f)df 



fPf{f)exp{-fl/AD)df = ^Ci/^^[fP^{f)], 



where Cs denotes the Laplace transform with respect to the variable s. We have the nor- 
malization condition Pi{l)dl = 1 and additionally 

If we now relate Pi{l) to g'i\l), 

ffl 2MV pn._ZD 

where Z = 2M/N = 2(i = 6 is the mean coordination number, we obtain the central 
theoretical result 

9i'\l) = {^Ci/AD[fPf{f)]- (7) 
D. Classification of Jammed Packings 

Jammed hard-sphere packings can be classified based on their density 0. However, such a 
classification is clearly not sufficient in order to distinguish between ordered and disordered 



15, hj. In 



(often called random, despite the shortcomings of such terminology) packings 
fact, packings can have various degrees of order in them, and for hard-sphere packings the 
dominant form of ordering is crystallization into variants of the FCC lattice. We can use 
a hypothetical scalar order-metric ip to measure the amount of order in a packing, such 
that ip = 1 corresponds to fully ordered (for example, the perfect FCC crystal), and = 
corresponds to perfectly disordered (Poisson distribution of sphere centers) packings. Very 
large jammed packings are thus classified based on their position in the density-disorder 
(0 — ijj) plane, as sketched in Fig. HI as taken from Ref . |l5| . A state of special interest is 
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Figure 2: A highly schematic plot of the subspace in the density-disorder {(j)—^) plane where strictly 
jammed three-dimensional packings exist. Point A corresponds to the lowest-density jammed 
packing, and it is intuitive to expect that a certain ordering will be needed to produce low-density 
jammed packings. Point B corresponds to the most dense jammed packing. Point MRJ represents 
the maximally random jammed state. This is the most disordered jammed packing in the given 
jammed category (locally, collectively or strictly jammed). We conjecture that the Lubachevsky- 



Stillinger packing algorithm 
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19| typically produces packings along the right (maximally dense) 



branch, and we do not know of an algorithm that produces packings along the left (minimally 
dense) branch. 

the MRJ state, representing the collection of maximally random jammed packings, believed 
to be closely related to the traditional but ill-defined concept of random close packing (RCP) 
in three dimensions if strict jamming is considered, and to have a density of about ~ 0.64 
in three dimensions^. Additionally, the perfect FCC crystal and variants thereof correspond 
to the most dense jammed packing, with ^ 0.74. This work will focus on these two points 
in the 4> — ip plane. However, it is possible to produce packings with intermediate amounts 
of order and densities, for example, by allowing partial crystallization. 



^ Contrary to popular belief, the traditional concept of RCP does not have a two-dimensional analog for 
monodispersc disks 0, ^3 ■ 
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III. COMPUTATIONAL RESULTS 

We use event-driven molecular dynamics js|] as the primary computational tool for our in- 
vestigations. This enables us to perform exact molecular dynamics on hard-particle packings 
very close to the jamming point, which is not possible with traditional time-driven molecular 
dynamics algorithms. The algorithm monitors a variety of properties during the computa- 
tional run, including the "instantaneous" pressure, as calculated from the total exchanged 
momentum in all interparticle collisions during a certain short time period At. By allowing 
the shape of the particles to change with time, for example, by having the sphere diameter 
grow (shrink) uniformly at a certain (possibly negative) expansion rate dD/dt = 27, one 
can change the packing density. If the change is sufficiently slow, the system will be in ap- 
proximate (metastable) equilibrium during the densification, and one can rather effectively 
gather quasi-equilibrium data as a function of density. 

Event-driven molecular dynamics in which the particles (quickly) grow in size in addition 
to their thermal motion at a certain expansion rate, starting from a random (Poisson) 
distribution of points, produces a jammed state with a diverging collision rate. This is the 
well-known Lubachevsky-Stillinger (LS) packing algorithm 18|,|l9|, which we have used and 
modified to generate all the disordered hard-sphere packings for this study. During the 
initial stages, the expansion has to be fast to suppress crystallization and maximize disorder 
jioj l. and delaying further discussion to later sections, we will assume that the disordered 
packings used in this study are representative of the MRJ state. It is important to note that 
the algorithm typically produces packings that have rattling particles, i.e., particles that 
do not have true contacts with particles in the jammed backbone of the packing, and can 
be removed without affecting the jamming category of the final packing. We will discuss 
procedures for identification of such rattlers in what follows. 

To our knowledge, no verification of the exactness of Eq. Q for disordered packings exists 
in the literature. The perfect FCC crystal is stable until rather low densities, and the pressure 
seems to be rather accurately predicted by the free-volume approximation in a wide range 
of densities around close packing. This has been observed in the literature and a suitable 



corrective term was determined 



20|. However, for disordered packings,^evious studies have 



identified a coefficient smaller than 3 in the numerator, namely 2.67 



2^. In Fig. El we 



numerically verify the validity of Eq. Q with very high accuracy for disordered packings. 
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Figure 3: The inverse of the "instantaneous" (averaged over several hundred colhsions per particle) 
pressure of a nearly jammed (isostatic) packing of 1000 particles, as it is slowly diluted (using 
a negative expansion rate for the particles in the molecular dynamics algorithm 7 = 10^^) from 
~ 0.627 until an unjamming particle rearrangement occurs. Up to this occurrence, the free- 
volume theoretical relation p = 5^^ is satisfied to very high accuracy. Rattlers have been removed 
from the packing. 

In Fig. m we show the change of the coefficient (the configurational heat capacity in units of 
A^A;) C = (1 — <f)/<f)j)p with density. Agreement with the theoretical C = d = 3 is observed 
sufficiently close to the jamming point, but with rapid lowering of the coefficient from 3 away 
from the jamming point. This is because for sufficiently large jamming gaps, contacts other 
than the M true contacts start contributing to the collisions, and the polytope-based picture 
we presented so far does not apply exactly. We demonstrate this in Fig. ID by showing the 
number of contacts which participate in collisions {active contacts) as the jamming point is 
approached. Our investigations indicate that previous studies did not examine at the range 
of densities appropriate for the theory presented above and did not properly account for the 
rattlers. 
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Figure 4: The coefficient C during a typical slow densification (expansion rate is 10^^) of a 10, 000- 
particle system, starting from an equilibrated liquid at = 0.5 up to jamming. The final packing 
has 259 rattlers, so the expected coefficient is 3 • 0.9741 ~ 2.92, a value which is shown with a red 
line. It is clear that close to the jamming point Eq. ^ is very accurate, but a marked lowering 
from a coefficient of 3 is seen for pressures lower than about 10®, likely explaining the coefficient 



2.67 reported in works of Speedy 



|. The inset shows the estimated "collisional" coordination, 
defined as the average number of different particles that a particle has collided with during a time 
interval of about 100 collisions per particle, during the same densification. The expected number 
6 • 0.9741 ~ 5.85 is shown (this number is not asymptotically reached exactly since some of the M 
contacts do not participate in collisions frequently enough to be registered during the time interval 
used), and we see that as many as 8 contacts per particle are active at sufficiently large jamming 
gaps. 

A. Disordered Packings 

We have verified in previous publications that LS packings are tvpically collectively 



jammed using a testing procedure based on linear programming 



Unfortunately, 



the linear programming library used in the implementation cannot really achieve the kind 
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Number of events (in 1000s per particle) 



Figure 5: The short-term ("instantaneous") pressure versus number of events (mostly binary col- 
hsions) processed by the molecular dynamics algorithm ^, corresponding to a total run of about 
half a million collisions per particle. For the 1000-particle packing the pressure is stable, but for 
the larger packings a systematic pressure leak is observed. 

of numerical accuracy that we require in this work, specifically that for packings which are 
jammed almost to within full numerical precision {6 = 10^^^ — 10^^^). Additionally, it can- 
not handle three-dimensional packings of more than about a thousand particles. Another 
test for jamming, which we have found to be reliable for the purposes of this work, is to 
take the final packing produced by the LS procedure and then run standard event-driven 
molecular dynamics on it for long periods of time (on the order of thousands to hundreds of 
thousands of collisions per particle) and monitor the "instantaneous" pressure. If the packing 
is jammed, this pressure will be stable at its initial value. However, if the packing is not 
truly jammed, we have observed that the pressure slowly decays with time, the slower the 
"pressure leak" the more "jammed" the initial packing is, as illustrated in Fig. Similar 



observations are made in Ref. [2lJ. In addition, we track the average particle displacement 
(from the initial configuration) and check to see if there is a systematic drift with time 
away from the initial configuration. The two tests always agreed: A pressure leak always 
corresponds to a systematic drift away from the initial configuration. 
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We have observed that LS packings densified to within numerical capabihty only pass 
this rigorous jamming test of having no pressure leak if during the final stages of the LS 
densification the expansion rate is very small compared to the average thermal velocity 
(maintained constant via a velocity rescaling thermostat j^) of the particles (about five 
orders of magnitude or less). Similar observations are made in Ref. j2l|. If the expansion 
rate is too fast, we have found that the packings jam in slightly hypostatic configurations, 
where there are not enough particle contacts to ensure jamming. In particular, some particles 
have 2 or 3 contacts (and of course rattlers are present) . In order for a set of balanced forces to 
exist (which as we discussed is a necessary condition for a packing to be locally maximally 
dense) when a particle has less than 4 contacts, these contacts must be in a degenerate 
geometric configuration, namely 3 coplanar or 2 collinear contacts. We have indeed verified 
that this is what happens in the hypostatic packings produced by the LS algorithm. The 
number of such geometric peculiarities increases with increasing expansion rate, and also for 
more ordered packings, as we discuss later. 

We illustrate the progress of the densification during the final stages of the algorithm in 
Fig. ini The figure shows, for several snapshots of the packing during the densification, the 
cumulative coordination number 



i.e., the average number of particles within a gap I from a given particle. We we will often use 
this quantity instead of (72(0- For the first time in the vast literature on random packings, a 
clear separation is seen between the delta-function contribution Z'^^^l), which becomes more 
localized around contact, and the background increase in the mean coordination from the 
isostatic contact value of Z = 6, which remains relatively unaffected by the densification. 
For small packings (A^ = 1000), the value of Z{1) is fixed at 6 for a remarkably wide range 
of gaps, as much as 9 orders of magnitude for the final packings. Fast densification is seen 
to lead to subisostatic packings in Fig. IHl leaving a certain fraction of the contacts "open". 
Stopping the expansion invariably leads to a decay of the macroscopic pressure for such 
subisostatic packings. 

By using heuristic strategies, we were able to find (slow) densification schemes which 
produced packings which are indeed ideally jammed within almost full numerical precision, 
at least for packings of = 1000 particles or less. In fact, the plateau in Z[l) was at 
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Figure 6: The cumulative coordination Z{1) (i.e., the integral of ^2(0) ^ ^ function of the gap 
tolerance /, for a sequence of snapshots of a 1000-particle packing during the final compression 
stages of the LS algorithm. For a sufficiently slow expansion (expansion rate is 10~^ times the 
average thermal velocity), the packing is clearly seen to jam in an isostatic configuration. A 
subisostatic configuration is found for fast expansion (expansion rate is comparable to the thermal 
velocity). The inset shows the properly normalized derivative of Z{1), right around contact, along 
with a comparison to our semi-theoretical prediction for g2\l), for a packing with 6 = 2.5 • 10^^^. 

exactly (up to a single contact!) an isostatic number of contacts, M = 3N — 2, for all 
the packings produced via a carefully guided LS algorithm. It is essential that here N is 
the number of particles in the jammed backbone of the packing jj], i.e., rattlers with 
fewer than 2 contacts have been removed from the packing. It seems that the algorithm 
produces packings with about 2.2% rattlers, and so the density of the disordered packings 
we look at is typically ^ 0.625 — 0.630, rather than the widely known ^ 0.64. Despite 
a concentrated effort and lots of expended CPU time, we have been unable to achieve true 
isostaticity for 10, 000-particle packings. This is illustrated in Fig. where it is clearly 
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seen that the pressure in the large packings does not remain constant over long periods of 
time (about a million collisions per particle). It is therefore not strictly justified to consider 
these packings within the framework of ideal jammed packings that we have adopted here. 
However, it is readily observed that over finite and not too long time intervals (for example 
several thousands of collisions per particle), the large packings conform to the predictions of 
the theory developed here. In particular, the coUisional forces form a balanced force network 
with essentially the same Pf{f) as the truly jammed smaller packings, and the pressure is 
given by Eq. (0)) with very high accuracy, where 6 can be determined, for example, via Eq. 
(jHl). We therefore believe it is justified to use the larger packings for certain analysis where 
better statistics are needed. 

The main goal of this work is to explore and explain Fig. IHl and in particular, to in- 
vestigate both the "delta-function", or contact, contribution g2^\ which should integrate to 
produce the isostatic average coordination Z = 2M/N = 6, and the "background" or near- 
contact for gaps from about 1006D to 10~^D. This latter one has already been observed 
in an experimental study of hard spheres |2,3], and in computational studies of stiff "soft" 



spheres |l 
and Ref. 



14\. These various studies find a nearly square- root divergence, g2 \i) ~ Vv^) 



24| observes that this is an integrable divergence and thus clearly separate from 



the delta function. Our results, shown in Fig. are the first unambiguous and precise 
separation of the two pieces of the pair correlation function around contact near jamming. 
Our numerical data has precision {6 < 10^^^) not previously attained, since such proximity 
to ideal jammed hard-sphere packings can only be achieved in a true hard-sphere algorithm, 
and at present only event-driven molecular dynamics seems to provide the required numeri- 
cal robustness. It is rather interesting that although graphs showing the hard-sphere (72(0 
the literature have clearly demonstrated a divergence in g2{i) near contact for at least three 



decades 



25j, this seems to never have been clearly documented or investigated. We are led 



to believe that researchers were under the false impression this divergence is a signature of 
the delta-function contribution, and thus expected it to further narrow and disappear at 
true jamming. 
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1. Delta- Function (Contact) Contribution 



A 







T 




1 



We first verify tliat our tlieory correctly predicts tlie sfiape of g2^\l)- In order to verify 
relation ((7j) numerically, a form for Pf{f) is needed. Force networks in p arti cle p ackings have 
been the subject of intense theoretical and experimental interest [3 |l^ ll^ls^ U^, and 
it has been established that Pf decays exponentially at large forces for a variety of models. 
The behavior of Pf for small forces has not been agreed upon, the central question being 
whether the infinite-system-limit Pf{0) is nonzero. No theoretical model has been offered 
yet that truly answers this question. Part of the difficulty is that the answer likely depends 
not only on the system in question, but also on the definition of /. In a true ideal collectively 
jammed isostatic packing, which is necessarily finite, all interparticle forces, must be strictly 
positive, and in fact are determined uniquely through Eq. (jS)), 



(8) 



without any mention of interparticle potentials or infiuence of external fields or loads like 
gravity, or thermal dynamics. The limiting probability distribution of these interparticle 
forces as the packing becomes larger, if it exists, can be positive at the origin, indicating 
that finite but large packings have limiting polytopes with a few extremely small faces, or 
equivalently, are very elongated along certain directions. We have numerically studied the 
form of Pfif) for almost jammed random packings of = 1000 and N = 10,000 spheres 
by using molecular dynamics to observe the collisional forces between first neighbors, and 
also by directly using Eq. (jHl) for the smaller packings^ (this offers better accuracy for small 
forces). The results are shown in Fig. [3 We clearly see a peak in P{f) for small forces, 
as observed in the literature for jammed packings of soft particles [2^, and it appears that 
there is a finite positive probability of observing zero interparticle force. We will return to 
this point later. 

The observed Pf{f) can be well fitted for medium and large forces by Pf{f) = [Af'^ + 
B)e^^'-^ , with a small correction needed to fit the small-force behavior, as used in Fig. El 
This small correction has a negligible impact on the Laplace transform of fPf{f), and in 



^ Efficiently inverting the rigidity matrix for very large three dimensional packings is a rather challenging 
numerical task which we have not yet tackled. 
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Figure 7: Computational data on the interparticle force distribution along with the best fit we 
could achieve. Packings of both 1,000 and 10,000 particles, using either molecular dynamics to 
average the coUisional forces, or inversion of the rigidity matrix, were used, consistently producing 
the same probability distribution. Comparison to other data in the granular-media literature is 
beyond the scope of this work. 

fact a very good approximation to g2\^) in Eq. (jT)) is provided by just using 

fiA R 

ifPfif)] = --T^ + --T^- (9) 

In the inset in Fig. we show a comparison between the g2\^) we observe computationally 
and the one given by Eqs. (|7j) and Q and the empirical fit to Pf{f) in Fig. [7| An essentially 
perfect agreement is observed. Our focus here is on small forces, however we do wish to note 
that our data cannot confidently rule out a Gaussian component to Pj for large forces and 
that a slight quadratic component does seem to be visible when Pf{f) is plotted on a log- log 
plot. 
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2. Near- Contact Contribution 



In Fig. ISlwe investigate the near-contact contribution to g2{l)- We have found that Z'^^\l) 
has a power law behavior over a surprisingly wide range of gaps, up to the first minimum 
of ^2 at / ^ 0.25L', Z^^\l) ^ ll{l/Df -^, as shown in the figure. Note that this range is too 
wide for 



240(1 + xy dx 

to be a perfect power law, where x = l/D, as used to fit numerical data in other studies 
(which have not investigated nearly as wide a range of gaps as we do here) Addition- 
ally, the observed exponent is clearly distinguishable from an inverse square root divergence 
in 5'2*''(0) as proposed in the literature 2^]. 

We do not have a theoretical explanation for this functional behavior of Z^^\l), however, 
the remarkable quality of the fit in Fig. |H1 hints at the possibility of a (simple) scaling 
argument. Some simple observations can be made by assuming that 

Z{x) = 2 + ax^-" for < X < /3, (10) 

where a is an exponent < « < 1, and (3 < 1 determines the extent of this power-law 
dependence. The corresponding pair correlation function of course exhibits an inverse power- 
divergence with exponent a, except when a = 1, when it is identically zero''. The exponent 
a clearly will depend on the amount of order present in the packing, i.e., the position of 
the packing in the density-order diagram of Fig. |2l We expect that it will increase with 
increasing order, since a ^ would indicate a constant (72(0 near contact, a signature of the 
ideal gas, while a 1 would indicate a clear distinction between the first and second shell 
of neighbors (i.e., a wide range of gaps with very few contacts) typical of crystal packings. 
Under the assumption that a power-law divergence in g2 is appropriate, an intermediate 
value of a between and 1, as we find numerically, is therefore expected. Some bounds 
on the range of possible a can be obtained from bounds on Z{x) derived from geometric 
constraints (for example, Z{x) < 13 for a certain range of x since the sphere kissing number 
is 12 in three-dimensions), but the exact value is not simple to predict^. 



^ Note that g^'^ (x) cannot have a simple-pole divergence since this would lead to a logarithmic divergence 

in Z^'''){x), which must be finite for all finite x. 
^ The three parameters a, f3 and a are thus not independent of one another. For example, requiring that 
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Figure 8: The near-contact Z^^\l) for a nearly-jammed 10, 000-particle packing, along with a power 
law fit for small gaps, shown in both a linear-linear scale and a log-log scale (inset). In this inset 
we also show a line with slope 0.5 (i.e., a square root dependence), which is clearly inconsistent 
with the numerical data. 

3. Away from Contact: Split Second Peak 

Although the primary focus of this work is on the behavior of g2{'r) around contact, it is 
instructive to also look at the split second peak of the pair correlation function, shown for 
a sample of packings of 10, 000 particles in Fig. El Only two clear discontinuities are seen, 
one at exactly r = -\/3-D, and one at r = 2D. The latter is very clearly asymmetrical, with 
a sharp decrease in g2 a.t r = 2D^. Although the first discontinuity is less pronounced and 
statistics are not good enough to unambiguously determine its shape, it appears that it also 
has the same shape as the second discontinuity, only of smaller magnitude. The split second 

92\^) > 1 ^{^) < 12 for < a; < /3 gives the weak constraints a(l — a) > 240/3^(/3 — 1)" and 
a{l3-lf~" <12- Z. 
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Figure 9: Computational data on the split second peak of g2{r) averaged over 5 packings of 10, 000 
particles. The values r = \/^D and r = 2D are highlighted, and match the two observed discon- 
tinuities. Also visible is the divergence near contact. The inset shows the probability distribution 
Pe{S) of bond-pair angles in the contact network of the packings, also revealing two divergences at 
9 = 7:/?, and at 9 = 27r/3. No peaks are observed at r = ^/2D or r = V^D, which are typical of 
crystal packings, indicating that there is no detectable crystal ordering in the packing. 



peak is of great importance because it is a clear signature of the strong local order in the 
first two coordination shells of the packing, and in fact observations have been made that 
along with the appearance of a peak in Pf{f) for small forces, the splitting of the second 
peak of g2 is a signature of jamming 2^. It is therefore important to try to understand the 
local geometrical patterns responsible for the occurrence of these structures in g2. 



4- Contact-Network Statistics 

The exact geometry of the jammed configuration Rj is determined (not necessarily 
uniquely) from its contact network, which as we have demonstrated is the network of first- 
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neighbor interactions and can easily be separated from further-neighbor interactions. Fig. 
[TUl shows the histogram of local coordination numbers as a function of the first-neighbor 
cutoff r, i.e., the histogram of the number of particles within distance (1 + t)D from a given 
particle. It is seen that for sufficiently small r (r < 10^^) the histograms are independent 
of the exact cutoff used (this is true down to r ^ 1005 or so, which can be as small as 
10^^^ in some of our packings). A number of particles having less than 2 contacts are seen, 
and these are clearly rattlers and we have removed them from consideration from all of the 
final packings we analyze here. We observe that such particles remain with fewer than 2 
contacts for a very wide range of r and are easy to identify. In some cases, however, we 
cannot unambiguously identify a handful of the particles as rattlers or non-rattlers. This is 
typical for packings which are not sufficiently close to their jamming point, packings which 
have been produced using fast expansion in the LS algorithm, or packings which are very 
large. It is safest to not remove such particles as rattlers. 

This work is the first time a clear look has been provided at the exact contact network of 
disordered hard-sphere packings. Previous studies have either used soft atoms, in which case 
the definition of a contact is not clear cut unless one carefully takes limits of a stiff interac- 
tion potential and therefore in such studies r has been typically set to correspond to the 
location of the first minimum in g2{f'), or have used Voronoi tessellations to define neighbors. 
Even studies which have actually used hard particles have resorted to such definitions un- 
suitable to investigating the jamming limit, mostly because the numerical precision required 
to separate the true contacts from the near contacts has not been achieved up to now 
Such investigations, the literature of which is too vast to cite, have found a plethora of local 



coordination patterns typical of polytetrahedral packing, inclu ding icosahedral order 



28|. 



We therefore attempted to do a similar shared-neighbor j28| analysis for the contact 
networks of our disordered packings, and look for local clusters reminiscent of polytetrahedral 
packing. Our procedure, based on looking at the contact network as an undirected graph, 
was as follows. For each particle, we extracted the subgraph corresponding to the first 
neighbor shell of the particle (this includes contacts between the neighbors), extracted its 
connected components, and counted the number of occurrences of a given subgraph (using 
graph algorithms that can test for graph isomorphism to form equivalence classes). The 
results were surprising. By far the most prominent patterns were a central particle contacting 
a chain of 1, 2, 3, 4 or 5 contacting particles. The chains were almost never closed, other 
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Figure 10: The probability distribution of local contact numbers as the cutoff used in defining 
neighbors is increased. Rattlers are clearly seen, and a relative maximum at Z = 6 is seen. Note 
that only one particle with 11 neighbors is observed, and very few have as many as 10 neighbors. 
No particle with 12 contacting neighbors has been observed in any of our packings, indicating a 
lack of crystallinity. 

than for chains of length 3 (which together with the original particle form a contacting 
tetrahedron), and this was itself rare. The probability of finding a chain of length n seems 
to decay exponentially, P{n) ~ exp(— 1.2n). This study found very few tetrahedra, and so 
polytetrahedral local ordering is certainly not apparent in the contact networks. We also 
performed the same analysis for a range of r's, all the way up to r = O.ID (which raises 
the average coordination significantly above 6), but still found the open linear chains to be 
the dominant pattern. We further attempted to include second neighbors in the analysis, 
however including all second neighbors led to very large subgraphs of a very broad variety, so 
classification was not possible. We further restricted our attention only to second neighbors 
which are very close to the given particle (within O.ID, for example), and this also found 
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very few tetrahedra. 

One of our goals was to determine if certain simple local coordination patterns are re- 
sponsible for each of the three features of (72 (^) we previously documented: the power-law 
divergence near contact, and the discontinuous, if not diverging peaks at r = and 
r = 2D. We had little success in accounting for the first one by restricting attention to only 
the first two neighbor shells in the true contact network. In particular, we looked at all the 
near contacts (for example, with gaps less than O.OID) and whether the almost contacting 
particles were in fact second neighbors in the contact network. Indeed most were, however 
the majority only shared one particle as a first neighbor, or two or three first neighbors which 
were not themselves first neighbors. It was therefore not possible to isolate one particular 
local geometry as responsible for the multitude of near contacts. An interesting quantity we 
measured is the probability distribution PeiO) of bond-pair angles 6 in the contact network, 
meaning the angles between two contact bonds of a given particle. This distribution is shown 
in the inset in Fig. 121 and shows divergences at 6 = 71/ 3 and 6 = 27r/3, which correspond 
to distances r = 2D sm{9/2) oi r = D and r = \/3D. Although there is no divergence at 
6 = IT, the corresponding distribution of distances does show a divergence at r = 2D. 

We had more success with a shared-neighbors analysis for the split second peak. This 
was because we could increase r and thus progressively relax the definition of first neighbor. 
We found that with increasing r, an increasing majority of particle pairs at a distance close 
to y/3D were second neighbors, and that an increasing majority of them shared two neigh- 
bors which were themselves neighbors. This corresponds to two edge-sharing approximately 
equilateral coplanar triangles, a configuration which has been suggested as being responsible 
for the first part of the split second peak j25j. Note however that we do not observe any 
discontinuity in (72 at r = 1.Q33D, which corresponds to two face-sharing tetrahedra, which 
is another configuration often mentioned in connection with the split second peak. A similar 
analysis for the peak at 2D indicated that the majority of particle pairs at this distance 
share one neighbor between them, which represents an approximately linear chain of three 
particles, a configuration which has long been known to be responsible for the second part 
of the split second peak of (72- 
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B. Ordered Packings 

In this work we have focused on disordered hard-sphere packings, and have found a 
multitude of unexpected singular features, such as a long power-law tail in g2\l), a nonzero 
Pf{f = 0), and a power-law divergence in g^\l)- It is important to realize that the properties 
we observe are not universal and will change as one changes the amount of ordering of the 
packings. In particular, dense ordered packings like the FCC crystal are not isostatic, and 
we have no theory that can predict the shape of . We therefore resort to a computational 
investigation of ordered and partially ordered sphere packings. 



1. Vacancy-Diluted FCC Crystal Packings 

It was the behavior of crystal packings around the jamming point that was the subject 



of Refs. |1| and |2| , and these works inspired this investigation. For crystal packings, there 
is no ambiguity in defining first neighbors, and the FCC packing has Z = 12 contacts per 
particle, which is twice the isostatic value. Therefore, the limiting polytope Px is not a 
simplex and, as argued in Ref. y, it is expected that for an FCC packing g2 \l) will have 
a single peak for small gaps. We indeed observe this computationally, for the first time, as 
shown in Fig. ITTl 

Furthermore, we have prepared vacancy-diluted FCC packings by removing a fraction p 
of the spheres from a perfect crystal, < p < 4 (here p = corresponds to the perfect 
crystal). The FCC lattice is composed of 4 interpenetrating cubic lattices. We obtain the 
vacancy-diluted crystal with the lowest density by removing one of these 4 cubic lattices 
(i.e., p = 1/4), as shown in the inset in Fig. ^2 This gives a packing with density of about 
(j) ^ 0.56 and mean coordination Z = 8 and is still collectively jammed. In fact, it is likely 



that more spheres can be removed with a more elaborate procedure la]. We can add back 
a randomly chosen fraction q = 1/4 — p of the previously removed quarter of the spheres, 
to obtain < p < 1/4. The delta-function contributions to g2 for several p's are shown in 
Fig. ^2 It is rather surprising that the pair correlation function for the p = 1/4 packing no 
longer shows a peak, but is monotonically decaying. In fact, by changing p one can obtain 
packings with g2 \0 that has zero slope at the origin. 

It is interesting to note that for the (vacancy-diluted) FCC packings g2\l) decays in a 



27 




Figure 11: The first-shell §2 (/) for a collection of FCC crystal packings with a fraction p of the 
spheres removed, starting with = 13, 500 particles. The inset shows the packing with most 
vacancies, where every 4**^ sphere is removed to form a cubic sublattice of vacancies (colored dark). 
Intermediate p's are achieved by randomly adding back some of the spheres to the sublattice. The 
density has been reduced by 5 = \/2 • 10^^^ from close packing. 

Gaussian manner, and in fact is perfectly fitted by a modified Gaussian, g2 \l) = (^^^ + 
Bl + C) exp [(/ — -D)^]. This fast decay is to be compared to the slow power-law decay for 
the disordered packings [c.f. Eq.®], hinting at possible connection to the stability of the 
crystal packings versus the metastability of the glass packings |8]. Additionally, we show the 
force distribution -P/(/) for these ordered packings in Fig. [121 illustrating that, in contrast 
with the disordered packings, very small forces are not observed. It would be interesting to 
know if the perfect FCC crystal can be vacancy-diluted to an isostatic packing and still be 
coUectivelly or strictly jammed, and what the corresponding force distribution would be. 
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Figure 12: The force probability distribution for the cohection of FCC crystal packings from Fig. 
111! For the pure crystal and the crystal with most vacancies, all of the particle pairs are identical 
and therefore the probability distribution would be a delta function if forces are averaged over an 
infinite time horizon. For the intermediate p's, multiple relatively broad peaks are observed. In 
contrast with the disordered case, very small forces are not observed. 

C. Partially Crystallized Packings 

As previously explained, the Lubachevsky-Stillinger algorithm can produce partially crys- 
talline sphere packings if a sufficiently small expansion rate is used and nucleation of crys- 
tallites occurs during the densification. This is demonstrated in Fig. ^[ where we show the 
evolution of the pressure during the densification of an initially liquid sample (i.e., a state 
on the stable equilibrium liquid branch), for a range of expansion rates 7. The slower the 
expansion is, the more crystalline the final packings become, as can be seen from the fact 
that the final density increases and from the evolution of the peaks in g2{f)i as shown in 
Fig. ^3 Additionally, the structure factor S'(k) shows more anisotropy and localized peaks. 

The packings shown in Fig. clearly have nucleated crystals, and so one may anticipate 
that there is a qualitative distinction between them and the "random" packings produced by 
suppressing crystallization. However, as demonstrated in Fig. slower densification leads 
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Figure 13: Compression of an initially liquid system with (p = 0.5 to jamming with several different 
expansion rates 7 (the mean thermal velocity is 1, in comparison). The pressure is plotted on a 
reciprocal scale (the tickmarks being equally spaced in equal increments of increasing in the 
usual direction), to highlight the expected linear relation Q near jamming. The pressure-density 



curves for the perfect FCC crystal [20|, the accepted fluid-solid coexistence region, and the widely 
known Carnahan-Starling equation of state for the fluid branch, are also shown for comparison. 
Sufficiently fast compression suppresses crystallization and leads to densities around 0.64 — 0.65, 
and slower compression allows for partial crystallization, typically occurring around (j) ^ 0.55, 
which is the end of the coexistence region (i.e., the density where the crystal necessarily becomes 
thermodynamically favored). This produces denser packings which exhibit more crystal ordering 
the denser they are. 

to larger densities and more ordered packings even if crystallization is suppressed and no 
visible nucleation occurs. This indicates that there is a continuum of packings from most 
disordered to perfectly ordered packings, so that one needs to be careful in interpreting 
results obtained from packings produced by just one, possibly non-trivially biased, algorithm. 
For example, Ref. j26| relates the occurrence of a peak in -P/(/) to jamming. However, as 
we show next, jammed packings do not necessarily exhibit this peak if they are sufficiently 
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Figure 14: The evolution of the peaks in g2{'r) as crystahine order is increased, for the packings 
from Fig. \VA\ The formation of peaks at distances typical of the FCC lattice, such as r = ^/2, 
is clearly seen. It is interesting to note that a peak is observed at y^ll/S ~ 1.91, which is a 
fifth-neighbor distance in the HCP (but not the FCC) lattice (a similar HCP peak at -^8/3 w 1.63 
is barely visible) 

ordered. 

For the sake of brevity, we will only briefly discuss some interesting features of g2 for 
the partially crystallized packings. Since the perfect FCC/HCP crystals have Z = 12, one 
expects that, as partial crystallization occurs, somehow the number of first neighbors per 
particle should increase from the isostatic value of Z = 6. However, this is not really so 
if one properly defines first neighbors via true contacts in the final jammed packing. In 
fact, if one plots Z{1) for partially crystallized packings (we omit this plot), a qualitatively 
similar curve to that shown in Fig. IHlis seen, with Z clearly close to the isostatic value of 
6. However, the background Z^^\l) shows a faster rise the more crystalline the packing is 
[consistent with a larger exponent a as defined in Eq. (fTUI) ]. so that indeed an increase of 
the cumulative coordination is seen for sufficiently large gaps. Additionally, we observe that 
nearly crystalline packings easily jam in noticeably hypostatic configurations, with a higher 
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Figure 15: Compression of an initially (metastable) liquid system with (p = 0.6 to jamming with 
several different expansion rates, as in Fig. lll^L For this range of expansion rates, crystallization 
is suppressed due to the large initial density and all final packings are apparently disordered and 
would be ordinarily identified as random, however, it is clear that slower compression leads to higher 
densities, and thus the final packings are not all identical, but rather, some are more ordered than 

n 

others, as can be verified by the slight increase in bond-orientational order metric Qg .15], for 
example. 

probability of observing particles with only 2 or 3 contacts and a less flat plateau in Z{1). 

All of these findings are readily explained. The basic premise, used widely in the granu- 
lar media literature, is that random perturbations to either the particle-size distribution or 
to the boundary conditions will break some of the contacts in an otherwise perfect crystal 
down to the isostatic value. This is because additional contacts in excess of Z = 6 imply 
special correlations between the positions of the particles, which one expects to destroy with 
random perturbations. Such random perturbations are provided in the case of partially 
crystallized packings by the fact that the crystallites need to jam against a partially amor- 
phous surroundings, and this induces complex strains that break some of the perfect-crystal 
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Figure 16: The evolution of Pf{f) as crystalline order is increased, showing the disappearance 
of the peak at small forces. The inset shows a log-log view of the plot, and is consistent with 
exponential decay for large forces. 

contacts.^ However, the geometric peculiarities of the underlying crystal remain; for exam- 
ple, there is a multitude of nearly collinear (in fact lines of aligned particles) or coplanar 
contacts, which leads to the occurrence of much more pronounced force chains (chains of 
large forces propagating along a nearly straight line) and a sharp increase in the probability 
of occurrence of small forces. We indeed observe this in Fig. El where we show that for 
sufficiently ordered packings there is no longer a peak in Pf{f) for small forces, but rather 
a monotonic decrease of Pf{f), apparently exponential for sufficiently large forces. This is 
in contrast to previous studies of the effect of order on force distributions in granular piles 
{29! , which did not register a significant impact of the ordering. However, these studies 
study the distribution of forces in granular piles and a direct comparison is beyond the scope 
of this work. 

^ We mention in passing that we have observed similar results by starting with a perfect FCC crystal, 
applying a small (but not too small) random strain, and then jamming the packings. This typically yields 
almost perfectly crystal packings which are nonetheless clearly frustrated by the random strain to have 
Z « 6. 
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IV. DISCUSSION 



The results presented in this work settle some long-standing questions and confusions in 
the literature. For the first time, we showed both theoretically and computationally how the 
delta-function portion of g2ir) is formed as jamming is approached, for a true hard-sphere 
packing. Our investigation focused on maximally disordered (MRJ) sphere packings with 
a packing fraction ?a 0.64 — 0.65. We presented the first true hard-sphere computational 
data on the power-law divergence in the near-contact portion of g2, in agreement with pre- 
vious observations in the literature for stiff soft spheres, but with a distinguishably different 
coefficient of —0.4. We confirmed that this divergence persists even in the true jamming 
limit for hard particles. We presented high-quality data on the probability distribution of 
interparticle forces -P/(/), especially focusing on small forces, demonstrating a maximum at 
small forces and a nonzero intercept at / = 0. A local analysis of the topology of the contact 
network found few traces of tetrahedra and an overwhelmingly complex local connectivity, 
and was successful in accounting for the structures responsible for the split second peak of 
g2{r)- A computational study of the delta-function contribution to g2{r) for vacancy-diluted 
FCC crystals showed a faster than exponential decay, unlike the slow power-law decay for 
the disordered isostatic packings. Finally, we investigated packings on the transition from 
maximally disordered to maximally ordered, and found that partially crystallized packings 
produced by our algorithm are still nearly isostatic despite having a higher density, and that 
Pfif) loses the peak for sufficiently ordered packings. 

This work has raised several important questions. The computational observations un- 
dermine the very applicability of the ideal jammed packing model to large (maximally) 
disordered packings of spheres, as produced by most algorithms in use today. First, a very 
unusual power law divergence in g2{l) is observed near contact, leading to a multitude of 
particle pairs just away from contact. Similarly, a power-law decay is seen in the contact 
part of (0- packings become larger, one can expect the tails of the two power laws 

to start overlapping by an observable number of contacts, blurring the distinction between 
true contacts and almost contacts. Even more troubling is the observation that there ap- 
pears to be a positive probability of observing a zero force in the contact network of the 
packings, indicating the presence of geometric degeneracies in the contact network. The 
above observations may explain why we have had trouble generating truly jammed packings 
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of = 10, 000 particles. However, we do not see a reason why very large but finite packings 
collectively jammed ideal packings could not be constructed. The question of what algorithm 
can produce disordered (and thus likely isostatic) packings which are jammed and devoid of 
some or all of the above peculiarities, as is the FCC crystal packing^^, for example, remains 
open. As usual, with each careful study the hard-sphere system provides more questions 
than originally posed or answered! 
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